#include "../includes/need-a-name.hh"
#include "includes/euler.hh"
#include <map>

Euler::Euler(){
  phynumber = 2+gridsetting::dimension;
  physical_variable_name.push_back("density");
  physical_variable_name.push_back("momentum1");
  if(gridsetting::dimension>=2)
    physical_variable_name.push_back("momentum2");
  if(gridsetting::dimension==3)
    physical_variable_name.push_back("momentum3");
  physical_variable_name.push_back("energy");  
  // TODO: transmit needs improvement
  transmit_equ_parameter();
  append_parameter("gamma",5.0/3.0);
}

Euler::~Euler(){
}

void Euler::get_equation_flux(int dim, double* dir, double* var, double* fequ){
  for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
    fequ[tmpi]=0.0;  
  if(dim==1)
    get_1d_fx(fequ,var);
  if(dim==2){
    if(std::abs(dir[0])==1.0)
      get_2d_fx(fequ,var);
    if(std::abs(dir[1])==1.0)
      get_2d_fy(fequ,var);
  }
  if(dim==3){
    if(std::abs(dir[0])==1.0)
      get_3d_fx(fequ,var);
    if(std::abs(dir[1])==1.0)
      get_3d_fy(fequ,var);
    if(std::abs(dir[2])==1.0)
      get_3d_fz(fequ,var);
  }
}

double Euler::get_pressure(double* var){
  int dim = gridsetting::dimension;
  double pressure;
  if(dim==1)
    pressure=(para["gamma"]-1.0)*(var[2]-1.0/2.0*std::pow(var[1],2)/var[0]);
  if(dim==2)
    pressure=(para["gamma"]-1.0)*(var[3]-1.0/2.0*std::pow(var[1],2)/var[0]-1.0/2.0*std::pow(var[2],2)/var[0]);
  if(dim==3)
    pressure=(para["gamma"]-1.0)*(var[4]-1.0/2.0*std::pow(var[1],2)/var[0]-1.0/2.0*std::pow(var[2],2)/var[0]-1.0/2.0*std::pow(var[3],2)/var[0]);
  return pressure;
}

void Euler::get_local_vmax(int index){
  int dim = gridsetting::dimension;
  double pressure;
  double var[memory::phynumber];
  for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
    var[tmpi]=memory::variable_tmp[tmpi][index];
  pressure=get_pressure(var);
  memory::local_vmax[index] = std::sqrt(para["gamma"]*pressure/var[0])*2.0;  
}

void Euler::get_1d_fx(double* fequ, double* var){
  double pressure;
  pressure=get_pressure(var);
  fequ[0]=var[1];
  fequ[1]=std::pow(var[1],2)/var[0]+pressure;
  fequ[2]=(var[2]+pressure)*var[1]/var[0];
}

void Euler::get_2d_fx(double* fequ, double* var){
  double pressure;
  pressure=get_pressure(var);  
  fequ[0]=var[1];
  fequ[1]=std::pow(var[1],2)/var[0]+pressure;
  fequ[2]=var[1]*var[2]/var[0];
  fequ[3]=(var[3]+pressure)*var[1]/var[0];      
}

void Euler::get_2d_fy(double* fequ, double* var){
  double pressure;
  pressure=get_pressure(var);    
  fequ[0]=var[2];
  fequ[1]=var[2]*var[1]/var[0];
  fequ[2]=std::pow(var[2],2)/var[0]+pressure;
  fequ[3]=(var[3]+pressure)*var[2]/var[0];  
}

void Euler::get_3d_fx(double* fequ, double* var){
  double pressure;
  pressure=get_pressure(var);    
  fequ[0]=var[1];
  fequ[1]=std::pow(var[1],2)/var[0]+pressure;
  fequ[2]=var[1]*var[2]/var[0];
  fequ[3]=var[1]*var[3]/var[0];    
  fequ[4]=(var[4]+pressure)*var[1]/var[0];
}

void Euler::get_3d_fy(double* fequ, double* var){
  double pressure;
  pressure=get_pressure(var);      
  fequ[0]=var[2];
  fequ[1]=var[2]*var[1]/var[0];
  fequ[2]=std::pow(var[2],2)/var[0]+pressure;
  fequ[3]=var[2]*var[3]/var[0];
  fequ[4]=(var[4]+pressure)*var[2]/var[0];    
}

void Euler::get_3d_fz(double* fequ, double* var){
  double pressure;
  pressure=get_pressure(var);        
  fequ[0]=var[3];
  fequ[1]=var[3]*var[1]/var[0];
  fequ[2]=var[3]*var[2]/var[0];
  fequ[3]=std::pow(var[3],2)/var[0]+pressure;
  fequ[4]=(var[4]+pressure)*var[3]/var[0];        
}
